{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d54447c6-cc17-4cf8-b47a-0936b8031709",
   "metadata": {},
   "source": [
    "## Script to produce Appendix Figure H1 in Carleton et al 2022: \n",
    "#### Constructing a hybrid Shared Socioeconomic Scenario (SSP) to approximate an income trajectory that is endogenous to climate change, following Burke, Hsiang, and Miguel (2015)\n",
    "\n",
    "Panel B shows the mean difference in income between SSP2 and a hybrid socioeconomic scenario in which SSP2 projected income is replaced by SSP3 projected income only for the poorest 60% of countries in 2010. \n",
    "\n",
    "Run the notebook sequentially to produce Figure H1. Charts are stored in `OUTPUT/figures/Figure_H1_SSP_compare`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "13bd8515-3fc3-4981-a535-562fdf23cff4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cf8fefb-e6ed-4fbe-9183-fbef8cd32292",
   "metadata": {},
   "source": [
    "#### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "id": "27fd3dde-dc4f-425f-827a-2cd515000992",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define paths\n",
    "REPO = os.getenv('REPO')\n",
    "DB = os.getenv('DB')\n",
    "OUTPUT = os.getenv('OUTPUT')\n",
    "\n",
    "# inputs\n",
    "gdppath = '{}/2_projection/2_econ_vars'.format(DB)\n",
    "poppath = '{}/3_valuation/inputs/interpolated_pop_exp'.format(DB)\n",
    "\n",
    "# output\n",
    "output_dir = '{}/figures/Figure_H1_SSP_compare'.format(OUTPUT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "50c13ec4-4fbe-4ba3-a365-a1e884559974",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define utility functions\n",
    "\n",
    "def wquantile(x,q):           \n",
    "    xsort = x.sort_values(x.columns[0])\n",
    "    xsort['index'] = range(len(x))\n",
    "    p = q * x[x.columns[1]].sum()\n",
    "    pop = float(xsort[xsort.columns[1]][xsort['index']==0])\n",
    "    i = 0\n",
    "    while pop < p:\n",
    "        pop = pop + float(xsort[xsort.columns[1]][xsort['index']==i+1])\n",
    "        i = i + 1\n",
    "    return xsort[xsort.columns[0]][xsort['index']==i]\n",
    "\n",
    "def reset_index(df):\n",
    "    '''Returns DataFrame with index as columns'''\n",
    "    index_df = df.index.to_frame(index=False)\n",
    "    df = df.reset_index(drop=True)\n",
    "    #  In merge is important the order in which you pass the dataframes\n",
    "    # if the index contains a Categorical. \n",
    "    # pd.merge(df, index_df, left_index=True, right_index=True) does not work\n",
    "    return pd.merge(index_df, df, left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ae58fff-6c3e-4e77-9243-587da43060cd",
   "metadata": {},
   "source": [
    "### set parameters and toggles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "58379919-741c-44b7-8ea6-de8cc9e8591f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" Creates dataframe and chart showing time series of difference in iso-level GDPpc of one\n",
    "    SSP scenario relative to another of 5 income quintiles divided based on a reference year\n",
    "\n",
    "    Distinct steps in this process include:\n",
    "        - identify income quintiles\n",
    "        - create timeseries of average GDPpc under each scenario\n",
    "        - take ratio of comparison to reference\n",
    "        - plot output\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "df1: SSP you are plotting difference of \n",
    "df2: reference SSP\n",
    "flatten: True if setting the top 2 quintiles to 0\n",
    "year: year which quintile breakdown is based; default is 2010\n",
    "logvar: False if using GDPpc in $ (default); True if using log GDPpc\n",
    "savefig: saves to OUTPUT/figures/Figure_H1_SSP_compare\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "SSPtest = 'SSP3'\n",
    "SSPref = 'SSP2'\n",
    "flatten=True \n",
    "year=2010\n",
    "logvar=False\n",
    "savefig=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "7bb2ac7d-ba23-4caf-b289-e5f0ef39f37d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import data \n",
    "df1l = pd.read_csv(os.path.join(gdppath, 'iso_gdppc_low_{}.csv'.format(SSPtest)))\n",
    "df1h = pd.read_csv(os.path.join(gdppath, 'iso_gdppc_high_{}.csv'.format(SSPtest)))\n",
    "\n",
    "df2l = pd.read_csv(os.path.join(gdppath, 'iso_gdppc_low_{}.csv'.format(SSPref)))\n",
    "df2h = pd.read_csv(os.path.join(gdppath, 'iso_gdppc_high_{}.csv'.format(SSPref)))\n",
    "\n",
    "dfpop1 = xr.open_dataset(os.path.join(poppath, 'interpolated_pop_exp_{}.nc4'.format(SSPtest)))\n",
    "dfpop2 = xr.open_dataset(os.path.join(poppath, 'interpolated_pop_exp_{}.nc4'.format(SSPref)))\n",
    "\n",
    "# extract only relevant pop data\n",
    "dfpop1 = dfpop1.squeeze('ssp').to_dataframe().reset_index()[['iso', 'region', 'year', 'pop']].groupby(['iso', 'year']).sum().reset_index()\n",
    "dfpop2 = dfpop2.squeeze('ssp').to_dataframe().reset_index()[['iso', 'region', 'year', 'pop']].groupby(['iso', 'year']).sum().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "15c5ccdc-3c68-465e-8ba3-1aa49ae9c5a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename gdp cols\n",
    "df1l.rename(columns={'0':'low_gdp'}, inplace=True)\n",
    "df1h.rename(columns={'0':'high_gdp'}, inplace=True)\n",
    "df2l.rename(columns={'0':'low_gdp'}, inplace=True)\n",
    "df2h.rename(columns={'0':'high_gdp'}, inplace=True)\n",
    "\n",
    "# combine low and high IAM in each GDP SSP file\n",
    "df1 = df1l[['iso', 'year', 'low_gdp']].merge(df1h[['iso', 'year', 'high_gdp']], left_on=['iso', 'year'], right_on=['iso', 'year'])\n",
    "df2 = df2l[['iso', 'year', 'low_gdp']].merge(df2h[['iso', 'year', 'high_gdp']], left_on=['iso', 'year'], right_on=['iso', 'year'])\n",
    "\n",
    "# create avg gdp per ssp\n",
    "df1['gdppc'] = (df1['low_gdp'] + df1['high_gdp'])/2\n",
    "df2['gdppc'] = (df2['low_gdp'] + df2['high_gdp'])/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "db7d5632-e121-43e4-90ea-42351e34ff93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# quintiles identified in df1 will equal df2 \n",
    "\n",
    "# take only 2010 and relevant columns\n",
    "df_temp = df1.loc[df1['year']==2010, ['iso',\"gdppc\"]].reset_index(drop=True)\n",
    "\n",
    "# drop duplicates (which are currently small countries just filled in by global average)\n",
    "df_temp = df_temp.drop_duplicates(subset=['gdppc']).reset_index(drop=True)\n",
    "\n",
    "# divide into 5 quintiles\n",
    "qlabs = [\"Poorest 20%\", \"20-40 Percentile\", \"Middle 20%\", \"60-80 Percentile\", \"Richest 20%\"]\n",
    "df_temp['qtile'] = pd.qcut(df_temp['gdppc'], 5, labels=qlabs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "29595fff-e128-498b-9554-68e17c53ae57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge quintiles into gdp dfs\n",
    "df1 = df1.merge(df_temp[['iso','qtile']], how='left', left_on='iso', right_on='iso')[['iso', 'year', 'gdppc', 'qtile']]\n",
    "df2 = df2.merge(df_temp[['iso','qtile']], how='left', left_on='iso', right_on='iso')[['iso', 'year', 'gdppc', 'qtile']]\n",
    "\n",
    "# merge pop data into gdp dfs\n",
    "df1 = df1.merge(dfpop1, how='left', left_on=['iso', 'year'], right_on=['iso', 'year'])\n",
    "df2 = df2.merge(dfpop2, how='left', left_on=['iso', 'year'], right_on=['iso', 'year'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "49f2fa94-bb33-4b13-8c54-74772d4fae54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate population weights by year, qtile \n",
    "df1['pweight'] = df1['pop']/df1.groupby(['year','qtile'])['pop'].transform('sum')\n",
    "df2['pweight'] = df2['pop']/df2.groupby(['year','qtile'])['pop'].transform('sum')\n",
    "\n",
    "# multiply out gdp by weight\n",
    "df1['wgdp'] = df1['gdppc'] * df1['pweight']\n",
    "df2['wgdp'] = df2['gdppc'] * df2['pweight']\n",
    "\n",
    "# take annual means by quintile\n",
    "df1 = reset_index(pd.pivot_table(df1, values='wgdp', index='year', columns='qtile', aggfunc='sum'))\n",
    "df2 = reset_index(pd.pivot_table(df2, values='wgdp', index='year', columns='qtile', aggfunc='sum'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "a38660fd-efee-456b-b961-cbaa6efcf97e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate ratio between df1/df2\n",
    "d = {}\n",
    "for q in qlabs:\n",
    "    d[q] = (df1[q] / df2[q] - 1)*100\n",
    "\n",
    "# set highest 2 quintiles to 0 if flatten == true\n",
    "if flatten == True:\n",
    "    d[\"60-80 Percentile\"] = d[\"60-80 Percentile\"] * 0\n",
    "    d[\"Richest 20%\"] = d[\"Richest 20%\"] * 0  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "id": "0e715d30-8f5d-4f90-8dee-f5857afd2457",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAAEWCAYAAABc752tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABfIklEQVR4nO2dZ3hVxdaA3xVC751QQm8BkkhHpV6KiKIIIk0Q9CKKXMTPXi5YQUEviCKCAoogooB0QbrSe+8QkN47pK7vx+yEAyQhgVMSMu/z7OecPXvPrLX3OWftOWvWrBFVxWKxWCy+w8/XClgsFktaxxpii8Vi8THWEFssFouPsYbYYrFYfIw1xBaLxeJjrCG2WCwWH2MNsY8RkY4iMtfXengTEVkkIs/dYd1AEbkkIuncrde9goj0F5GXfa1HWkNE/iMiA+6k7j1jiEWknYisFJHLInLCef+iiIhzfIyIRIjIRWfb4nxhc7q08YyIRDs/9AsiskFEHvGk3qo6TlWb3kldEcklIqNE5JhzTbtE5A2X448513BBRE6JyHwRKeEc6ycikc61nhORZSJSxznWUEQ2O+WnRWSKiBRxywUn/xrDRKRx7L6qHlTVbKoa7WY5z4jI33dR31OfRQsR+dspPyYiI0UkeyJ65Ac6A9/e6bXc6zj3+6ebykRE9onItnjOXyQiKiIhN5X/7pQ3cIpGAJ1EpEBydbonDLGI/B8wBBgIFAIKAj2AB4AMLqd+pqrZgfxAV6A2sFREsrqcs1xVswG5gO+BiSKSx+MXcWf8D8gGVARyAi2BvQAiUgb4Efg/51hJYBgQ41L/F+da8wN/A5OdB9c2oJmq5gIKA7uBb5KikIj43/VVpU489VnkBD7CfA4VgaKY73lCPAPMUtWr7rowb+LD7089oABQSkRqxHN8F+YBB4CI5MXYj5OxZap6DZjtel6SUdVUvWG+qJeB1rc5bwzw0U1l2YGjwEvO/jPA3y7HswIKVI+nvWeApcBQ4DywA/iXy/E8wGjgCHAW+D0BvW6WqZiHyG6n3teAJFB3C/B4AsfaABsSuR/9gJ9c9is5svPddF5GoD+wLZG2woA3gE1AOODvfEmXAeeAjUADl/MXAc8570sDC4DTwClgHJDLOTYWY6yuApeA14ESjp7+QDtgzU269AGmueg+CDgIHAeGA5nj0b8icA2IduScc/lu/Yj5sR0A3gX8fPVZOMeeADYn0tYCoJPLfm5ghnMNZ533RZ1jt7t/eYHpwAVgNeaB8Hcisn8FjmF+D0uASi7HMgOfO/fxPOZhk9nl83zW+ZyWYDqI7zrnnnA+g5xOO5mAn5zvyzlHr4Iuv6V9wEVgP9AxKffbKRvlfPcmA1/ddGwR8F/gEJDOKXsJ0zk5xI3f7Y7AwsRsUXzbvdAjroP5wU1NbkVVvQj8CdS9+ZjzZH4O88PcnUATtTAffD6gL6YXE9t7HgtkwfyoCmB6TEnlEaAGEAK0BZolcN4K4GMR6SoiZW86tg6oICL/c1wN2RISJiIZMV/iQ6p6yikLFJFzGCP4KvDZbXRuD7TA/JMoCMzE/HDzOPUnOX+bbxGPMfSxPb5imB8Kqvo05sf5qBp3xM06TAPK33TtHYDxzvtPgXJAKFAGKIL5Qd2Aqm7HPPyWO3JyOYeGYoxxKaA+pqfTNYHr99hncRP1gK0J1QeqADtd9v0wHYLiQCDm8/zKOXa7+/c1ppNTCOjibIkxGyiL+b6vwxi2WAYB1YD7Md+J17nxH0F9zOffDHP9zwANMfc+m4vOXTCfSTHMg6IHcNX5V/sl0FzNv977gQ230RcAEcmCeViOc7Z2IpLhptOOYP4pxroRO2MeEDezHfO7TR7JtdwpbQM6AcduKovtiV0F6jllY7ipR+yUDwD+dHmiRjl1T2F+XI0TkPuM8+GIS9kq4GkgAPMly50E/Z/h1h7xgy77E4E3E6ibGXgbWAtEAnswX8TY47Wd+icxPb4xQDaXXkGEc60nMD2pavHIyIPp7dZO5BrCgG4u+28AY286Zw7QxaWH8VwCbT0OrL+p7cYu+yWce+Tv7P8E/Nd5XxbTG8qCMfCXgdIudesA+5P4OaTD9O6DXMqeBxb58LNogunVlkvks4gEKiRyPBQ467Kf0P1L57RV3uXcRHvEN8nJ5XxOOTEPg6tASDznxX6epVzK5gMvuuyXd3TxB7phft/BN7WT1bl/rYnnX89N5/bjxn8gnZzPxR/TqTsHtHI5vgjTKesE/Ozos8s5dnOPuCwQnZR75LrdCz3i00A+V9+Sqt6vpldzmtv7wYsAZ1z2V6hqLlXNp6q1VXVeInUPq3P3HQ5genbFgDOqejY5F+LCMZf3VzA9gltQ1auq+omqVsP0DiYCv8b2ylV1haq2VdX8mF5/PeAdlyYmOtdaQFUbqeraeGScAX4Apt7Gf/ePy/viwJPOANM5p2f9IOYBdQMiUkBEJojIYRG5gDEM+RKRczPjMb1xML2531X1CsbXmgVY66LDH055UsiHGV844FJ2APN9uQVPfxYiUtu51jaquisRvc9iXG6x9bKIyLcicsC5v0uAXHI96iSx++fPjZ+r6/sbEJF0IjJARPY6csKcQ/mcLROOzzwBXNsuzK333R/zT2ss5qE+QUSOiMhnIpJeVS8DT2F6yEdFZKaIVEhEnitdMPc/SlXDMe6J+Hr/k4FGQC9Hj/jIjnG9JIt7wRAvx/RcHktuRecvYmPgrzuUXcQZUIklENNL/gfIIyK57rDdZKOqF4BPMD2DkvEcX435IlW+g+b9MX83cySmgsv7fzA94lwuW1ZVjS+0p79TN1hVc2B6Ha73VOOp48pczIM4FGNQYv9Wn8L0wiq56JBTzYDY7fSPrR+JeajEEggcvo0+bv8sROQ+jBuhm6rOv83pmzDumFj+D9ODq+Xc33qxzTqvCd2/k5h/h0Vd2iqWiNwOmN9gY0wvuISLnFOYfwGlE6nvev+PcOt9jwKOq2qkqr6vqkEY98MjOINjqjpHVZtgHvg7gJGJyDPKiRTFGNdOTlTKMYyb4mERuaFD4DygZgMvkLAhrogZE0kWqd4Qq+o54H1gmIi0EZFsIuLnfLGyxldHRDKKSDXgd0wPYvQdii8A/EdE0ovIk5gPYZaqHsV8YMNEJLdzvF6iLd0BIvKeiNQQkQwikgnojflbtVNEHhSRf8eG0ji9g5YYd8vt2n1CRMo79zE/8AXGXXDmdnUdfgIeFZFmTk8pk4g0cL70N5MdZ4BMTIjcazcdP47xE8aLqkYBv2EiCfJgfP6oagzmh/g/l3tQREQS8rcfB4rG+gbVhMdNxPh9s4tIceAV59puwYOfRWVMT76Xqk6/3fnALIy/NZbsmAfSOad33tf15ETuXzTmYdHP6VVXIPFogOyYDtFpzD+RT1xkxGAGw74QkcLOd6KO4w+Pj5+BPiJS0uksfYKJKolyfOxVnB79BczDMlpECopIS8dXHI75TiUlxPFpTEREeYzbJhTzIDvE9X8KrrwN1FfVsATaq4/57SeLVG+IAdQM4ryCGQA4gflRfYvxVS5zOfV1EbmIcUX8iPHn3e/8rbkTVmJ8QqeAjzF/G087x57GfEl2ODq9fIcyEkMxD5FTmF5EE6CFql7CGIGWwGYRuYT5MU/h9oNuYP5+/4HxF27G+LtbJVkp1X8wvaO3MT2rfzAGNr7v2/tAVczfuZmYH78r/YF3HffCqwmIHI/pif3qGJZY3sD4alc4f5fnYX5w8bEAMwh2TERiB8l6YfzM+zCj/OMxBiU+PPVZ/B/GTfC9mDjjSyKS2GDdj5jeXGZnfzDGfx075vFHPHUSun8vYXq3xzA9wJ8xRi4huQcw/xi2cetD5lXMd2k15vf3KQnbn1GOvCWY6IdrmM8CzMDhbxgjvB1YjHk4+mHu1RGn/frAiwm070oXYJiqHnPdMBE2t7gnVPWIqsYbb+48gB/GuPKShdzo4rQkFRF5BjPg9KCvdbFYXBGRT4ATqjrYze1+ChRS1dtFT6RJRKQXUExVX09u3bQafG+x3LOo6tvuaMdxR2TA9GRrYGJ972hqelpAVYfead0U6ZoQkWIislBEtovIVhHp7ZT3c0bXNzjbw77W1WK5h8mOcRVdxvjLP+cO4vUttydFuiZEJAAIUNV1YubVr8XEl7YFLqnqIF/qZ7FYLO4kRbomnKiDo877iyKynQTiNy0WiyW1kyJ7xK6IyVC1BBNz+QpmBtQFYA3wf/FNmhCR7kB3gKxZs1arUCGpcd0Wi8WSNNauXXvKmaBz16RoQ+zEEC4GPlbVySJSEBOGo8CHGPdFt8TaqF69uq5Zs8bzylosljSFiKxV1eruaCtFDtYBiEh6YBIwTlUnA6jqcVWNdgnWr+lLHS0Wi8UdpEhD7Ewb/h7YrqpfuJS75ipohUk9aLFYLKmaFDlYh0no/jRmJtIGp+xtoL0zdVkxSUWe94VyFovF4k5SpCF2phBKPIdmeVsXi8Vi8TQp0jVhsVgsaQlriC0Wi8XHJMkQO2n8ujrv84vILTlWLRaLxXJn3NYQi0hfTDrBt5yi9CSQk9VisVgsyScpPeJWmFyql8Hk48RlKRaLxWKx3B1JMcQRzrpsCuBkwLdYLBaLm0iKIZ4oIt9iFhz8N2aVg9uuBWWxWCyWpHHbOGJVHSQiTTCJdspjlt7+0+OaWSwWSxrhtobYiZD4K9b4ikhmESmRyOJ5FovFYkkGSXFN/IpZPDKWaKfMYrFYLG4gKYbYX1UjYnec9xk8p5LFYrGkLZJiiE+KSMvYHRF5DJMT2GKxWCxuIClJf3oA40TkK0winn+Azh7VymKxWNIQSYma2AvUdlbLEFW96Hm1LBaLJe2QlKiJjEBroATgb3K2g6p+4FHNLBaLJY2QFNfEVOA8Zkn7cM+qY7FYLGmPpBjioqr6kMc1sVgsljRKUqImlolIFY9rYrFYLGmUpPSIHwSeEZH9GNeEAKqqwR7VzGKxWNIISTHEzT2uhcVisaRhbuuaUNUDQDGgkfP+SlLqWSwWiyVp2BU6LBaLxcfYFTosFovFx9gVOiwWi8XH2BU6LBaLxcckGjUhZj7zL0AF7AodFovF4hESNcSqqiLyu6pWA6zxtVgsFg+QFNfEChGp4XFNLBaLJY2SlAkdDYEeIhKGiZywM+ssFovFjdiZdRaLxeJj7Mw6i8Vi8TF2Zp3FYrH4GDuzzmKxWHyMnVlnsVgsPsbOrLNYLBYfk6AhdhYNRVUHAb8Bk7g+s26od9SLV6+HRGSniOwRkTd9pYfFYrG4i8TC15YDVUVkrKo+TQqYWSci6YCvgSbAIWC1iExT1W2+1cxisVjunMQMcQYR6QLcLyJP3HxQVSd7Tq0EqQnsUdV9ACIyAXgMSNAQR4cf4cL+fi4lcv1VQBAQp0zkxnNEAD/zKn5mcz3fYrFY3EBihrgH0BHIBTx60zEFfGGIiwD/uOwfAmrdfJKIdAe6A4RWDoBorttf1HmrzsnX62lcA65t3aqEqmnFbKDq8h5x9i0WiyVpJGiIVfVv4G8RWaOq33tRp8SIryt6i81T1RHACIDq1atrjjL9AIi8coWLhw5xafduru7YwbW9e4ncv5+oI0fQ48eRM2fIGBlJZiALkAfIlhHIAdElcuFXuyJSvhAEZAG/i3DxEFw6DFdO3KqGfyYo1hAqdIAyj0EGG/FnsdxbvOO2lpIyxXmCiLwLBKpqdxEpC5RX1Rlu0yLpHMLM8oulKHAkqZXTZ8lCnnLlyFOuHLRocctxVSX83Dku/vMPFw4eZNOaNZycPJnMmzcTePIcgWtXkiUmxpxbrBjSsCE07AUPPQh5MxijHLud2wd7fofZT4N/Zij1KFRoDyWbg3/Gu7sLFovlnkJMiHAiJ4j8AqwFOqtqZRHJDCxX1VAv6HezLv7ALuBfwGFgNdBBVbcmVKd69eq6Zs2au5J7+dgx9s2axb4ZM7gwZw4BV64QKEJxf38yRkYCoCVLOobZ2YoUAY2BI8th+zjY9StcPQUZc0KZJ4xRDmwIfkl5FloslpSGiKxV1epuaSsJhniNqlYXkfWqep9TtlFVQ9yhQHIRkYeBwUA6YJSqfpzY+e4wxK5ER0Rw6K+/2D97NvtmzsRvxw6KAaWzZKFIVBTpIyLMieXKQbdu8OKLkD07REfCwfmw42fYMwUiLkKWAlCuLVRoB4XrOIOBFoslNeBtQ7wM0wNdqqpVRaQ08LOq1nSHAp7G3Yb4Zs7t30/YnDnsnz2bf+bNI9eVKwT6+VExRw4KnjuH5s2LvPIKvPQS5MhhKkVehf2zjFHePxOirkH2QGOQy7eDAqE2MsNiSeF42xA3Ad4FgoC5wAPAM6q6yB0KeBpPG2JXosLDObJsGfv/+IPdkyeTac8e/pU/PwEnT0Lu3PDyy/Cf/0CuXNcrhV+AvVONUT7wJ8REQe7yULEDVHkOshX2iu4WiyV5eM0Qi4gf0AaYD9TGRC2sUNVT7hDuDbxpiF2Jjoxk/VdfsaxvX/JdvUrTkiXJt3u36RX36mWMcr58N1a6ehp2TzZG+Z9F4JcOyrSC0J5QtJ7tJVssKQhv94iXqGo9dwjzBb4yxLFcPnaMxa+/zraxYyleoADNSpQg++rVSObM0KMHvPoqBATcWvHcXtg4HLZ8D9fOQt5KEPoiVOwEGXN4/0IsFssNeNsQvwdcxazmfDm2XFXPuEMBT+NrQxzL4WXLWNCrF8fXrSOoWjUaFipE5j/+AH9/6NoVXn8dSpa8tWLkVdg5ATYMg+NrIH02COoEIS9AfrtalcXiK7xtiPfHU6yqWsodCnialGKIAWKio9kyahR/vf02V0+fptZTT1E7fXrS//ILREdDhw7w5psQFBR/A8dWG4O8c4IZ4Cv8AIS+AGXb2Nhki8XLeNUQp3ZSkiGO5drZsyx7/33Wf/UVGbJlo17v3lQ5dw6/776DK1fg8cfhnXegegKf8dUzsHUMbPwGzu2BzPmgcjcIfh5ypYrno8WS6rGGOBmkREMcy6lt21j0yiuEzZlDnvLladS3LyV27IChQ+HsWXjuORgwAPLmjb8BjYED841B3jsNNBpKNIPgHlD6ETtZxGLxIO40xHYGgQ/JFxRE69mzaTVjBqrKbx068NuqVZyeMwdeew1Gj4YKFcyrM7X6BsQPSjSBxybDvw/A/e/DqS0wrRWMLAkrP3HyYFgslpSMNcQ+RkQo3aIFz2zeTIMvvuDo8uWMqVOHhZGRRK1cCeXLmxl69erBpk0JN5S9CNT5L/w7DFpOgTwV4O93YEQxmPU0HFlh0sZZLJYUR5JcEyLSEogNYVusqtM9qpUbScmuifi4cvIkS997j43ffktArVq0nDiR7PPnmx7yuXNmQki/ftdn6SXG6R2w4WvY9oOZUl2gqolJrtAO0mfx9KVYLPc0XnVNiEh/oDcm+fo24D9OmcUDZMmfnybDh9Pyt984tXUrY2vU4J/SpWHnTnj2WRg82PSSf/rp9j3cvBXgX0Ph+cPwr2EQHQ5zn4URRWHR/8HZ3V65JovFkjhJCV/bBISqaoyznw5Yr6qpIog1tfWIXTm1bRtTW7Xi3N69NBg0iKq9eyNr1kDPnrB6NdStawb2QpKYf0kVDv8F67+GPZPNdOriTSDkRTu4Z7EkE18M1uVyeZ/THYIttydfUBCdVq2i9KOPsrBPH2Z26EBExYqwYgWMHAnbt0PVqsYwn0nC/BoRM1X60V/g3wfh/g/g9Pbrg3vLP4BLSU7vbLFY3ERSesTtgQHAQkyuiXrA26r6s+fVu3tSc484Fo2JYdVnn/H3O++Qp2JFHps0iTzly5sQt7594euvTSKhjz6C7t0hXbqkNx4TBXtnmBC4A3NB0pkVRSp3M6FwtpdsscSL1+OIRSQAqIExxCtV9Zg7hHuDe8EQx3Jg3jxmtG9PdHg4zUaNonybNubA5s1mEG/RIggOhiFDoEGD5As4txc2fmsmi1w9CVkDIKgzVO4Kecq78UosltSPtwfr5qvqUVWdpqpTVfWYiMx3h3BL8ijeuDFPr1tH3qAgpj/5JAtfeYXoyEioUgUWLIDffoPz580KIR06wIlkxhDnKg31P4PnD5kQuILVYc0gGF0Bxt8Pm0ZA+HnPXJzFkoZJ0BCLSCYRyQPkE5HcIpLH2UoANkmuj8hRrBjtlizhvl69WPu///FLgwZcPHTI+H9btzZ+4759YdIkqFgRxoxJfvxwugxQ9nFoNc0Y5XoDIeI8/Pk8DC8EMztC2J8QE+2JS7RY0hwJuiZEpDfwMsboHub6CsoXgJGq+pU3FLxb7iXXxM3smDCBOf/+N/6ZMtFi3DhKNG16/eD27fDvf8PSpdCoEQwfDmXL3rkwVZP9bcsY2DEews9BtqJQqTMEdYE85e72ciyWVIW3s6/1UtWh7hDmC+5lQwxwZudOprVpw6mtW6n97rvc37cvfrGDdTExJrrijTfg2jWTSOj11yHjXWZqi7oGe6aaiSJhc0zOi4A6UKkLlH8KMuW66+uyWFI6NulPMrjXDTFA5JUrzH/pJbaMHk2xBg1oMX482VyTzR89Cn36wC+/mNwV33xzZ4N58XHpCGz7yRjl09sgXUYo29os81Ssvl0Q1XLPYg1xMkgLhjiWLT/8wLwXXiBD9uw8+uuvFKt308Iqs2ebmOP9+6FLFxg4EPLnd49wVTi+1kRcbB9nXBc5S0GVZ43rInsR98ixWFII1hAng7RkiAFObd3KtNatOb9/Pw+NGUPF9u1vPOHKFfj4Y2OEs2UzaTafew783NhzjbxqZu5t/s6svSd+ULK5iU0u9YgZDLRYUjneDl8TEekkIv919gNFpKY7hFvcT75KlWi/bBkBtWszs0MHVn76KTc8bLNkMYZ440YzNfr556FOHVi3zn1KpM8MFTtC24XQbTfUfBNOrIdpreHborDoVePGsFgsQNKmOA8D6gCxXauLwNce08hy12TOk4c2c+dSoV07/nrzTf7s0cPEG7tSsaKJPf7pJzhwAGrUMG6Ls2fdq0zuMvDgxyZfcqsZUORBWD8ExlSC8XVg00gIv+BemRZLKiMphriWqvYErgGo6lnA/rdM4fhnzEiLceOo9dZbbBoxgimPPEL4+ZsmY4hAx46wYwe89JIJcStXziSid7fLys8fSrUwSeyfPwz1P4eIC/BndxObPLuzcWNoPAnwLZZ7nKQY4kgn45oCiEh+wP5aUgHi50fdTz6h2XffcXDBAn5+8EHOHzhw64m5cplp0WvXGkPcrZuJqti+3TOKZSkA1V+BLlugwwozjXrPVJjYEL4vCyv7w+VUM4veYrlrkmKIvwSmAAVE5GPgb+ATj2plcStVnn2W1rNnc+HgQcbVqsXRlSvjPzE0FP76C777zuSvCAmBd9+Fq1c9o5gIBNSCJsOhx1FoPhZyBMLfb5uVRaa1vh6nbLHcwyQ16U8F4F+Y2XXzVdVDXSX3k9aiJhLj1LZtTHnkES4fPcpDP/xAhbZtEz755El49VX48UcoVQq++gqaN/eOomd2moiLLaPh2mnIUcKEwVXqasPgLCkGb0dNDAHyqOrXqvpVajLClhvJFxREx5UrKVi9OjOeeoplH3xAgg/i/Pnhhx/MgF769PDww9CmDRw65HlF85SH+gONL7nFBMhVCpa+ByMDYUpL2DvdpO+0WO4RkuKaWAe8KyJ7RGSgiLjlCWDxDVny5+fJefMI6tyZZX37MqN9eyKvXEm4QsOGJtTto49g5kwzM2/QILg5CsMT+GeECk/Bk/Oh2y6o/hocWwW/t4SRJYxxPr/f83pYLB4myRM6nExsrYF2QKCq3kUGGe9hXRPxo6qsHjiQJW++ScGqVXl86lSyF7nN3/79+6FXL2OQK1c2Celvnr3naaIjYd9047rY/wegENjYTKku87gx3haLF/DFUkkAZYAKQAlghzuEW3yHiFDz9ddpNXUqZ3bu5Kfq1TmyYkXilUqWhOnT4fff4eJFqF8fnn7a5LLwFunSQ9kn4IlZJja5Tj84txtmtoNvC8PCl+HkZu/pY7G4gaRkX/sUeALYC0wEJqvqOc+r5h5sj/j2nNyyhd9btuTS4cM0GTGCyl263L7SlSvQvz989pnJ5vb++yYWOX16zyt8MxoDB+abXvLe3yE6AgrVNAN85dtBxhze18lyz+PtNJg9gN9U9ZQ7BHoba4iTxtXTp5neti0HFyygWp8+1B848Ho6zcTYvdss0/THHxAUZFaVbtTI8wonxJVTsP0n2PI9nNoC/pmh3JPGKBepa0LmLBY34BVDLCIVVHWHiFSN77iqujE5wQ1yBwKPAhGYXnhXVT3nrAyyHdjpnLpCVXvcrj1riJNOTFQUi/7v/1j35ZeUfvRRWowfT4Zs2W5fUdW4LF5+2fiR27aFL76A2/mcPYkqHFttDPKOnyHiIuQqY9bfs9ngLG7AW4Z4hKp2F5GF8RxWVfVIt0dEmgILVDXKcYugqm84hniGqlZOTnvWECefDd98w/yXXiJ/SAitpk+//SBeLNeuGVdF//7g7w/9+pnesi/cFa5EXoZdk2DLKDi02GSDK9HMxCWXbmkH+Cx3hLddE5lU9drtyjyBiLQC2qhqR2uIvcu+2bOZ3rYtGXPmpNW0aRSsGu8fowQq7zMGeOZMqFTJTAZxVyL6u+XsHpPEfssYuHQIMuWBCh1MT7nAfdZ1YUky3o6aWJbEMk/QDZjtsl9SRNaLyGIRqZtQJRHpLiJrRGTNyZMnPa/lPUip5s3psHQpki4dP9ety67Jk5NRudT16IrLl6+vKn3kiMf0TTK5y8ADH8K/w6D1HCjeFDaPhJ+qwdhQWD3IrDpisXiRxFwThYAiwE9AB64vHpoDGK6qFe5YqMg8oFA8h95R1anOOe8A1YEnVFVFJCOQTVVPi0g14HegkqommkPR9ojvjsvHjvF7q1YcXbGCup98Qs0330SS02u8etUkn//0U+Oi+O9/oXdvyJCCEvhdOws7JsC2H+HoCuO6CGwMQU9D2VaQPquvNbSkQLzlI+4CPIMxhq6W7CIwRlWT0UVKplJGdg/gX6oa77QvEVkEvKqqiVpZa4jvnqhr15jz7LNsHz+eih070nTkSNJnzpy8RvbuNevmTZ8O5cubbG/NmnlG4bvhzC7YPha2jYULB4wRLvsEVOwEgf8CvyREkljSBN72EbdW1UnuEJYUROQh4AugvqqedCnPD5xR1WgRKQX8BVRR1TOJtWcNsXtQVVb278/f77xDoRo1eGzKlKQP4rkya5bpEe/ZAy1bmuiK0qXdr/DdojFw+G9jkHdONLmTswYYf3JQJ8gfYv3JaRyvr1knIi2ASkCm2DJV/cAdCsQjaw+QETjtFK1Q1R4i0hr4AIgCooG+qjr9du1ZQ+xe9kydysxOnciQPTuPT5lCQK1ayW8kPNz0iD/8ECIi4JVX4J13zBp6KZGoa7Bvhlmtev8siImEvJXMclAVO0CO4r7W0OIDvN0jHg5kARoC3wFtgFWq+qw7FPA01hC7n5ObN/P7Y49x6fBhGg8fTpWuXe+soSNH4M03YexYCAgwvuROndy7kKm7uXoadv1qjPKRpaasyIPGKJd7EjLn9a1+Fq/hbUO8SVWDXV6zYaY5N3WHAp7GGmLPcPX0aaY/9RQH58/nvl69aPD556S703jhFSuMu2LVKqhVC4YNg+SEy/mK8/th+3jYPg7ObDfLQQU2hvJPmQREmXL5WkOLB/F2+Frs8gxXRKQwEAmUdIdwS+olc968tPnjD6r/3/+xfuhQfm3cmMvHj99ZY7Vrw/LlMGYMhIWZhUx79YJz59yosQfIWRJqvwPPbIVO66Da/8GZHTCnKwwvaHInb/vJLo5quS1J6RG/BwzFrNDxNWbtuu9U9T3Pq3f32B6x59k+fjxznnuOTLlz03LSJArXrn3njZ07B++9Z3rF+fObmXpPP516BsZip1bvnAi7JsLFfyBdRijxEJRvC6UfhQzZfa2lxQ14fbDORXBGIJOqnr/tySkEa4i9w4mNG5naqhUXDx2i8ddfE/zvf99dg+vWQc+exm3xwAMm93FIiHuU9RYaA0dWGJ/yrl/h0mFjlEs2N/5ka5RTNd6KI34isYqejCN2J9YQe4+rZ84ws317wubOJeSFF2g0eDDp7mbiRkyMcVe88QacOQMvvggffAC5c7tNZ6+hMXB4mTHIu38zs/die8rl2picFzZdZ6rCW4Z4dCL1VFW7uUMBT2MNsXeJiY7mr7feYvXAgRStW5dHf/uNrAUK3F2jZ84Yd8Xw4ZAnj4mu6No1ZUdXJEa8RjmDmTBS5nEo/RhkLehrLS23wWeuidSINcS+Yfv48cx59lky58tHy0mTCKhZ8+4b3bDBJJ9fuhSqV4cvv4Q6de6+XV+iMXBkOeyeDHumOGvwCRS+3xjlMo9B7lSxKlmaw9urOBcUke9FZLazHyQiqSKG2OI7KnboQPtly5B06ZhQty6bv//+7hsNDYW//oKffjIxyPffD507p4xkQneK+EGRB6DB5/DsXui8Eer0hagrsOQ1GFUOxlSCv96Go6uM4bbccyQlamI2MBqTkCdERPyB9apaxRsK3i22R+xbrp4+zYx27Tgwbx7B3bvT6Msv8c/ohvy/ly7BJ5/A55+bZELvvGNyWWTKdPu6qYXzYbB3KuyZCoeWgEabadalW5qecrFGNpeyD/H2hI7VqlpDRNar6n1O2QZVDXWHAp7GGmLfExMdzd/vvMOqTz+lYPXqtPztN3IWd9O04L174dVXTcrNkiWNYX788dQT7pZUrp6B/TONUQ77wyS7T5/NJLgv8xiUaA5Z8vlayzSFtyd0XBaRvJj4YUSkNpBqwtcsvscvXTrqDRjAY1OmcHbXLsZWrcr+OXPc03jp0jBlCvz5J2TJAk88AY0bw5Yt7mk/pZA5j0nL2fI3ePEUtJppplUfWQazO8M3BWB8HVjxERxfZ10YqYyk9IirYiZ0VAa2APkxq2Zs8rx6d4/tEacszu7ezdTWrTm1ZQt13nuPOv/9b9IWKU0KUVEmsuK//4ULF6BHD7O6dN57OP+DxsDxtbBvlukxH1ttyrMWMr3kks2heBM73doDeM01ISLpgP9gDHF5THL4naoa6Q7h3sAa4pRH5JUrzOvZk61jxhD4r3/RYtw4shZ0Y7jW6dPGGA8fDjlymLXzXnzR92vneYPLxyFsjskSFzYHws+BpIPCdYxRLtEcCoTee64bH+BtH/EiVW3gDmG+wBrilMvmUaOY37MnGXPlosX48QQ2bOheAVu2mAG8efNMMvpBg6BFi7RjhGKi4OhKY5T3zYKTG0x51kJmIklcbzkVTpBJAXjbEH8M5AR+AS7HlqvqOnco4GmsIU7ZnNy0ielt23J2927q/Pe/1H73Xfe5KsDkfpgxwwzo7dpl/Meffw7Bwe6TkVq4dNTpLc+GA3Od3rIfBNQ2g34lmkHB6nYVkiTibUO8MJ5iVdVG7lDA01hDnPKJuHSJeS++yLaxYynWsCEtxo0jW0CAm4VEGFdFv35w/ryZmffBB1C4sHvlpBZiokxccths2P+H8TOjZlXrwMZQoqlZWDVHMV9rmmKxM+uSgTXEqYctY8Yw78UXyZA9Oy3GjaN448buF3LmDHz8MQwdanzGr71messpdXUQb3HlFBycZ3rMB+ZeX8k6T8XrRrlYfbuQqgvWECcDa4hTF6e2bmV627ac3r6d2u++y/3//S9+/v7uF7R3L7z1Fvz6KxQqZKIrunUDT8hKbajC6W2OUf4TDi2GqKvgl95MvS7R1PiWC1RN024Ma4iTgTXEqY+Iy5eZ/9JLbB0zhqJ16/LwuHHkKOahv8grVpge8dKlULGi6S3fixNC7oaoa2Yh1QN/Qtjc64N+mXKbREXFmxh3Rq5SPlXT21hDnAysIU69bPvpJ/584QXSZcjAQ6NGUeaxxzwjSNXMzHvrLdi5E2rWNBne3B3Fca9w+TgcnG8M84E/TZ5lMCuWBDY2xjmwEWTJ71s9PYy3B+t6AuNU9Zyznxtor6rD3KGAp7GGOHVzdvduZrRrx/F16wh54QUaDBpE+ixZPCMsKgp++MEM6B06BE2amB5yjRqekXcvoApnd8GBeY4bYxGEOxNv8wc7RvlfULTePZcE39uG+Ja8Eq55J1I61hCnfqLCw/n7nXdY8/nn5A0K4pGffya/J8PPrl6Fb74xSYVOn4ZWrUyEReXKnpN5rxATZaZYH5xnes2Hl0J0uFlYtVDN673lgDqpPmGR11dxBkLUOdGZbbdJVSu5QwFPYw3xvUPY3LnM7tKFa2fOULd/f6q9/DLiyeTwFy7A4MEm7vjiRWjf3vSWy9r8wEkm8qrJh3FwPvyzwEzB1hjwzwSFHzBGuWgDKFTdJMdPRXjbEA8ESgDDMYl/egD/qOr/uUMBT2MN8b3FlZMnmfPcc+ydNo3ARo1o/sMPZC9a1LNCT5+GgQNNIvqICJMD+b33TLY3S/IIP29Seh5cYAzzSSdljX8WE5FRrAEUrQ+FaqT4HrO3DbEf8DxmFWcB5mJWcY52hwKexhriew9VZfP337Pw5ZfxS5+eJt98Q4V27Twv+NgxM4g3fDhER5twt3fegcBAz8u+V7lyCg4vgX8Wme3UZlPun8m4L4rWN/HLhWpB+sy+0zMevB41ISIZgIpADCbpT4Q7hHsDa4jvXc7u2cOsp5/m6IoVlG/blsbDhpHZG5nWDh+G/v1hxAiz/9xzJuLCUyF2aYmrp+HQX2bQ75/FcHIjoMZtUaimMcxF65necwbfTsLxdo+4BcYtsRfTIy4JPK+qs92hgKexhvjeJiYqitUDB7K0b18y581L0xEjKP3oo94RfvCgMcjff2/ijp95xszUK1PGO/LTAtfOmRjmQ4uNS+P4WrNSiaSDglWhyIPG11zkAZPMyIt42xDvAB5R1T3OfmlgpqpWcIcCnsYa4rTBiY0bmfX005zavJmgTp1oOGQImfPk8Y7wAweMy2L0aIiMhCefhDfegPtSRWBR6iLikhn8O7TEbMdWmagMMHHMsUa58AOQN8ijM/+8bYiXqGo9l30BFruWpWSsIU47REdEsOLjj1n5ySdkypuXJt98Q9lWrbynwNGjMGQIDBtmoiyaNoXXX4dGjexMPU8RHWHC5Y4sNQb68FK4ctwcy5DDZJYrfL/ZAmpBxhxuE+1tQ/wNUByYiImaeBLYCSwFUNXJ7lDEU1hDnPY4sWEDf3TtyokNGyj35JP8a+hQ9yaevx3nzpkBvcGD4fhxqFbNuCxat7a5LDyNKpzfB0eWG8N8ZCmc2uIsHSWQr7JJkh9Qx7zmLmtSgd4B3jbEoxM5rKrazR2KeApriNMm0ZGRrBk0iGXvv0/6LFlo8MUXVOrSBfFmz/TaNRg71iSk37XLhLv16WNScKb1bG/eJPyCcWEcWWYM9NHl12f/ZcwJBauZgcBCNaBgDcheNEn/YGyuiWRgDXHa5vSOHcx97jkOL11KYKNGNPn2W3J7ezAtJgamTTOxyMuWQe7cZj29l15Ku/mQfYnGwJkdcGQFHF9tJpmc3GhmBYIZ9CtY3THMzms8eTOsIU4G1hBbNCaGjSNGsOSNN4iJiKD2e+9R49VXSZfBBzO5li83M/UmTzZuivbtTS85NNT7uliuE3XNGONjq69vZ3bgLF4P2QONQY41zgWrIpnzWEOcVKwhtsRy6cgR5vfqxe7Jk8kbFESTb7+l6IMP+kaZvXvNwN6oUXD5MjRoYBY4ffzxtLHIaWog4qIZCDy+5rpxPr8v7rC8ijXEScUaYsvN7J0xg3k9e3Lx4EEqd+tGvU8/JUu+fL5R5uxZMzFk2DATl1yokJkg8txzULy4b3SyJMzVM3BiHRxbjdR+222GOMHhQhGpJSIbReSSiCwXkSB3CLwdItJPRA6LyAZne9jl2FsiskdEdopIM2/oY7n3KP3II3Tdto0ar7/Oth9/ZFT58mwcMQKNifG+Mrlzm5jjffvMIqfVq5vUmyVLmhWnp00z6TktKYPMeaB4Y6j1llubTbBHLCJrgLeAJUBL4DlV9bjxE5F+wCVVHXRTeRDwM1ATKAzMA8rdLueF7RFbEuPU1q3Me/FFDi1ZQqGaNWk0ZAiFa9f2rVIHDsB335kZe0ePQpEi8OyzZrN5LVIM7hysSyyAzk9V/1TVcFX9FfB1uv3HgAmOPvuBPRijbLHcMfkqVeKpRYt4eOxYLv7zD+Pr1GHW009z8fBh3ylVvDh8+KExyJMnQ5UqZr9ECXj4YZgyxczgs9wzJGaIc4nIE7FbPPue5CUR2SQio5wVQQCKAP+4nHPIKbsFEekuImtEZM3Jkyc9rKoltSMiBHXqxLO7dlHr7bfZ+euvfF+uHMvef5+Iy5d9p1j69CYp/ezZxnXxzjuwcSM88YRJMPTGGyY+2ZLqScw14bGJHCIyD4gvQ8c7wArgFCZu5EMgQFW7icjXwHJV/clp43tglqpOSkyWdU1Yksu5/ftZ8sYb7Pr1V7IVLsyDn3xCpaef9mwS+qQSFQV//AEjR8LMmSYd54MPmoRDbdpAzpy+1jDNkGbiiEWkBDBDVSuLyFsAqtrfOTYH6KeqyxNrwxpiy51yeOlSFvbpw7HVqykQGkr9gQMp3rixr9W6ztGj8OOPJgRu1y7ImBFatoSOHaF5c/BFnHQawls+YkSkvIh8LiIznW2QiJRzh+BEZAa47LYCtjjvpwHtRCSjiJQEygKrPKmLJW1T5IEH6LhiBS3GjePa2bP82qQJvzVvzslNm3ytmiEgwLgnduyAVavg+edh0SITi1yokNlfvNjM7LOkaBILX6sDLAIuASOAkcBlYJGIeHJY+TMR2eysldcQ6AOgqlsxiYe2AX8APVPLKiGW1Iv4+VGxQwe67dhB/YEDObpiBT+EhjL7mWe4cPCgr9UziJiVpocMMUnrZ80yg3rjxpmJIoGB8OqrsG6dSYpjSXEk5iOeDXyqqotuKq8PvKmqzT2v3t1jXRMWd3L1zBlWDRjAui+/BCD0xRep9fbbvpsQkhiXL8P06TB+vBnwi4qCcuWgXTsztbpCqkgpnmLxio9YRHaparxuCBHZqarl3aGAp7GG2OIJLhw8yNK+fdn244+kz5qV6v/3f1Tr04eMOdyX79atnD5tQuF+/tm4L1QhJMQY5bZtoVQpX2uY6vCWIV6rqtUSOLZOVau6QwFPYw2xxZOc2raNpe+9x+7Jk8mUOzfV+vThvl69yJQrl69VS5gjR+C332DCBJOECIxRbtXKbFWq2ET2ScBbhvgEMCG+Q0BbVfVipu07xxpiizc4tmYNyz/8kL3TppExZ07u+89/qPbyy95brulOCQszPeUpU2DpUtNTLlXKDPg9/jjcfz+k89xyQ6kZb0VNvAasjWdbA7zuDuFpmXTp0hEaGkrlypV58sknuXLlilflh4WFMX78+HiPbdiwgTp16lCpUiWCg4P55Zdf4o7t37+fWrVqUbZsWZ566ikiIsyC3pMmTaJSpUrUrVuX06dPA7B3717aeWOZ+xRAoerVaTV1Kk+vW0dgo0as+PBDRhQvzuI33uDy8eO+Vi9hSpSAV16Bv/4y4XAjRkD58vDVV1CvnonM6NYNpk4FL39H0xSqek9v1apV05RI1qxZ49536NBBP//887tqLzIyMlnnL1y4UFu0aBHvsZ07d+quXbtUVfXw4cNaqFAhPXv2rKqqPvnkk/rzzz+rqurzzz+vw4YNU1XVOnXq6IULF3TEiBH65Zdfqqpqu3bt4tpJa5zYvFmnt2+vg/z89H+ZMumfL76o5/bv97VaSef8edUJE1Tbt1fNmVMVVDNnVm3ZUvW771SPHfO1hj4HWKNuslMJLqAlIg8CpVT1R2f/NyD2f9ZHqrrA848Jz7Pg5Zc5sWGDW9ssEBpKo8GDk3x+3bp12bRpE2fOnKFbt27s27ePLFmyMGLECIKDgxMs79evH0eOHCEsLIx8+fIxZMgQevTowUEnrGrw4ME88MADLF68mN69ewNmOu+SJUt488032b59O6GhoXTp0oU+ffrE6VOu3PUx2sKFC1OgQAFOnjxJzpw5WbBgQVxPukuXLvTr148XXngBPz8/wsPDuXLlChkzZuSvv/4iICCAsmXLuuGOpj7yV67MI+PH88D777Pqs8/YNHIkG7/9lort21Pj9dfJX6WKr1VMnBw54KmnzBYRYeKRp00zPeNp04wPuVYtM4HkkUegcmXrV74LEnNNvI9xQ8RSHuOu6Id1TbiNqKgoZs+eTZUqVejbty/33XcfmzZt4pNPPqFz584ACZYDrF27lqlTpzJ+/Hh69+5Nnz59WL16NZMmTeK5554DYNCgQXz99dds2LCBv/76i8yZMzNgwADq1q3Lhg0bbjDCN7Nq1SoiIiIoXbo0p0+fJleuXPg7C2AWLVqUw05ynL59+9KsWTPmzZtH+/bt+eijj3jvvfc8ddtSDbnLlqXZyJH8e98+qv7nP+yeMoUfgoOZ1KIF/yxZgqaGuN4MGaBJExg61CQiWr8e3n/fhMO9/TYEB5tERT16GEN96ZKvNU59JNRVBlbftD/Z5f1Sd3XJPb2lVNeEn5+fhoSEaEhIiL700ksaHh6uoaGhunfv3rhzihYtqufOnUuwvG/fvtqvX7+48vz588e1GRISooULF9YLFy5o//79tWbNmjpkyBD9559/VDVx10QsR44c0XLlyuny5ctVVfXEiRNaunTpuOMHDx7UypUr31JvzJgxOnjwYF2+fLm2bt1an3vuOb18+fKd3ah7jCunT+uyDz/Ur/Ln14GgY2vW1B0TJ2p0VJSvVbszDh82roonnlDNls24MDJkUP3Xv1QHDVLdulU1JsbXWnoE3OiaSMwQ707k2B53KeDpLaUaYlcfcSwhISG3GNzz588nWN63b18dOHBgXHnevHn1ypUr8crbtGmTDhgwQIsUKaLbt2+/rSE+f/683nfffTpx4sS4spiYGM2bN2+cP3rZsmXatGnTG+pdvnxZGzZsqBEREdqgQQM9f/68fvPNNzpixIjb3JG0RcTly7r+m290ZJkyOhB0RMmSumrQIL165oyvVbtzwsNVFyxQffVV1UqVjHkB1cBA1e7dVSdNUj13ztdaug13GuLEXBM7RKTFzYUi8giw0739cgtAvXr1GDduHACLFi0iX7585MiRI8Hym2natClfffVV3P4Gx/e9d+9eqlSpwhtvvEH16tXZsWMH2bNn5+LFi/HqERERQatWrejcuTNPPvlkXLmI0LBhQ3777TcAfvjhBx577LEb6n722Wf07t2b9OnTc/XqVUQEPz8/r0eFpHTSZ8lCaI8edNuxg5aTJpG9aFEWv/oqw4sWZe7zz6ecfBbJIUMGaNjQrFa9ZYtxY4wYAdWqmYkkrVtD3rxQty589BGsWWPzYMSSkIUGymAM7migl7ONAXZhVsbweW83KVtq6hGfPn1aW7ZsqVWqVNFatWrpxo0bEy2/uUd88uRJbdu2rVapUkUrVqyozz//vKqqvvTSS1qpUiUNDg7Wdu3a6bVr1zQiIkIbNWqkwcHB+sUXX9ygx9ixY9Xf3/8GN8f69etVVXXv3r1ao0YNLV26tLZp00avXbsWV+/w4cM39LInTpyoQUFBev/99+uJEyfcc+PuYY6vX69/PPus/i9TJh0IOv7BB3X7zz9rVHi4r1W7eyIiVBcvVn3rLdWqVa/3lvPmVW3dWvXrr1W3b09Vbgzc2CNONA2miGQEOgKVnKKtwHhVveapB4O7sRM6LKmNq6dPs2X0aDZ88w3n9+0jS4ECVHnuOYK7dyfnvbKg6IkT8OefZluwAP5x1nwICDC96titVKkUG42RZvIRuwNriC2pFY2JYf+cOWz85hv2zZyJqlKyeXNCnn+eUg8/jJ9/gtGnqQtVswLJ/PmwcKHZYifBBAaaDHING5rXEiV8qOiNWEOcDKwhttwLXDh4kE0jR7L5+++5fPQo2YsWpXK3blTu1u3e6SXHogrbt5vkRAsXmtdTp8yxWMNcv77ZfNhjtoY4GVhDbLmXiI6MZN+MGWz89lvC5s4FoESzZlTu2pUyLVvinymTjzX0ADExsG2bMciLFpnJJbGGuXBhs1TUAw+YLSQEvPRPwRriZGANseVe5XxYGJtHjWLLqFFcOnyYjLlyUeGppwjq3JnCdeogKdS3eteoGsO8eLFJVPT33xCbpD9rVqhd2xjl++837z20jp9XDbGIPICZTVcc8MdkX1NVTRUJTK0httzrxERHc3DBArb9+CO7Jk0i6upVcpUuTcVOnQjq1IncZcr4WkXPc/AgLFtmDPPSpWa165gY47aoXBnq1DGG+f77oUwZt7gzvG2Id2CWK1oLxC1NpKqn3aGAp7GG2JKWiLh4kV2TJ7Nt7FgOLlgAqhSqWZOgjh0p/9RTZC2YKrLX3j0XL5p1/JYuNTmXly+H8+fNsXz5TE+5dm0T41ytGuTPn2wR7jTEt41vA1a6K1bOF1tKjCM+ePCgNmjQQCtUqKBBQUE6ePDguGOnT5/Wxo0ba5kyZbRx48Z65jYzrQYOHKiAnjx5Mq7sk08+0dKlS2u5cuX0jz/+iLde/fr1tVy5chocHKz333+/7tixwz0Xlwz279+v48aNi9tfvXq19urVS1VVR48erT179vS6TvcSF/75R1d+9pmOCQnRgaCD/Pz016ZNdfOYMXrt/Hlfq+ddoqNVN29WHTFCtWtX1QoVrscyg2qxYqqPP676wQeqM2eqHj162ybxxhTnuBNgADAQqANUjd3cpYCnt5RoiI8cOaJr165VVdULFy5o2bJldevWraqq+tprr2n//v1VVbV///76+uuvJ9jOwYMHtWnTphoYGBhniLdu3arBwcF67do13bdvn5YqVUqj4sljUL9+fV29erWqqn777bf66KOPJkn3mJgYjY6OTvrFJkJi06ytIXYvJzZv1iVvv60jSpbUgaBfZMyov7dqpdt/+UUj0moekHPnVBcuVB04ULVdO9Vy5W40zgEBqi1aqL73nuqUKaphYTdMOHGnIU7K8GIt59W1C65AI7d0yX3Nyy+Dm9NgEhoKiaTBDAgIICAgAIDs2bNTsWJFDh8+TFBQEFOnTmXRokWASTPZoEEDPv3003jb6dOnD5999tkN04ynTp1Ku3btyJgxIyVLlqRMmTKsWrWKOnXqJKhPvXr1GOzoO3DgQCZOnEh4eDitWrXi/fffJywsjObNm9OwYUOWL1/O77//zi+//MLYsWPx8/OjefPmDBgwgL1799KzZ09OnjxJlixZGDlyJBUqVOCZZ54hR44crFmzhmPHjvHZZ5/Rpk2bW1Jx3nfffQwaNIgZM2bcoN/JkyfjTe9pSTr5K1cm/8cf8+BHH3F05Up2TJjAzokT2T1lCv5ZslDq4Ycp16YNpVq0IEO2bL5W1zvkzGlC4Ro0uF524YLJLrd+vVn1eu1as/Bq7FTsPHnM77uqe1eKu60hVtWGbpVouYGwsDDWr19PrVrmeXf8+PE4Ix0QEMCJEyfirTdt2jSKFClCSEjIDeWHDx+mdu3acfuuqSoTYvr06VSpUoW5c+eye/duVq1aharSsmVLlixZQmBgIDt37mT06NEMGzaM2bNn8/vvv7Ny5UqyZMnCmTNnAOjevTvDhw+nbNmyrFy5khdffJEFC0za6qNHj/L333+zY8cOWrZsSZs2bRgwYMANhjf2AXQzsek9H3zwQQ4ePEizZs3Yvn37be6sJT5EhMK1a1O4dm0afP45h5YsYeevv7J78mR2/fYb6TJmpESzZpR94glKP/poyl/qyd3kyHE9RjmWK1dg82ZjmGON9NChbhWbWGL4Tqr6k4i8Et9xVf3CrZr4imQkcHc3ly5donXr1gwePDjeJD4JceXKFT7++GPmOnGkrph/TDeSUBhTx44dyZw5MyVKlGDo0KEMGTKEuXPnct9998Xpt3v3bgIDAylevHicgZ83bx5du3YlS5YsAOTJk4dLly6xbNmyG5IEhYeHx71//PHH8fPzIygoiOPJXDpo3rx5bNu2LW7/woULXLx4kezZsyerHcuN+KVLR2DDhgQ2bMi/hg7lyLJl7PrtN3ZPmcLeadOQdOko1qABZVu1oszjj5O9SBFfq+wbsmQxSfBr1bpeFhlpkhy5icR6xFmdV/tt9wCRkZG0bt2ajh078sQTT8SVFyxYkKNHjxIQEMDRo0cpUKAAAF27dmX9+vUULlyYTz/9lP3798f1hg8dOkTVqlVZtWoVRYsW5Z/YefvOscKFC8erw7hx46he/brHSVV56623eP755284LywsjKxZs95w3s3GPSYmhly5csVlfLuZjBkz3lA/OcTExLB8+XIyZ86crHqWpOOXLh1F69alaN26NBw8mONr17Jr0iT2TJnC/JdeYv5LL1GoRg3KPPYYZR5/nLxBQfdunHJSSJ/erc0lmAZTVb91Xt+Pb3OrFmkMVeXZZ5+lYsWKvPLKjX84WrZsyQ8//ADcmGZy9OjRbNiwgVmzZlGlShVOnDhBWFgYYWFhFC1alHXr1lGoUCFatmzJhAkTCA8PZ//+/ezevZuaNWsmSa9mzZoxatQoLjkrLBw+fDhe10jTpk0ZNWpUXGrLM2fOkCNHDkqWLMmvv/4ad40bN25MVF5iqThvlhdfek+LZxARClWvTr3+/em2Ywddt2+n7iefIH5+/P3uu4ypXJnvy5ZlYZ8+HFy4kOjISF+rnOpJLB+xxUMsXbqUsWPHsmDBAkJDQwkNDWXWrFkAvPnmm/z555+ULVuWP//8kzfffDNZbVeqVIm2bdsSFBTEQw89xNdff026JC6H3rRpUzp06ECdOnWoUqUKbdq0iddQPvTQQ7Rs2ZLq1asTGhrKoEGDANPD/v777wkJCaFSpUpMnTo1UXnBwcH4+/sTEhLC//73vwTP+/LLL1mzZg3BwcEEBQUxfPjwJF2PxT3krVCBWm+9RccVK+hx+DCNv/mG3OXLs+Gbb5jYqBFf58vH1DZt2DxqFJeOHPG1uqkSO8XZYrHcERGXLnHgzz/ZP3s2+2bN4pIzKJw/JISSDz1EyebNKVynDunc6EtNSdhcE8nAGmKLxfOoKqc2b2bf7Nnsnz2bI0uXEhMVRfps2Qhs1IgSzZpRslkzcpUu7WtV3Ya3pzgXBD4BCqtqcxEJAuqo6vfuUMDTWENssXif8AsXOLhgAWFz5rD/jz+4EBYGQM5SpSjRtCnFmzQhsFEjMuXK5VM97wZvG+LZmOWS3lHVEBHxB9arahV3KOBprCG2WHyLqnJuzx7C5s4lbM4cDi5cSOSlS4ifHwWrVyewYUOKNWhAkQcfTFWTSbxtiFerag0RWa+q9zllG1Q11B0KeBpriC2WlEV0ZCRHV64kbO5cDs6fz7FVq4iJikLSpaNQ9eoUrV+fYvXrU+SBB8jooRSW7sDbhngR0Br4U1Wrikht4FNVrZ9oxRSCNcQWS8om4vJljixbxj+LFnFo8WKOrlpFTGQk4udHgdBQitarR9F69Sjy4INkuYMsaZ7C24a4KjAUqAxsAfIDbVQ1Vaz3bQ2xxZK6iLxyhaMrVvDP4sUcWrKEI8uXE+3M0sxToQJF69alSN26FHnwQXKWKOGziSXuNMS3jSNW1XVAfeB+4HmgUmoxwikZEeHpp5+O24+KiiJ//vw88sgjgMklMWDAgHjrZkvAj/bMM8/w22+/AdCgQQOS8wB67bXXqFChAsHBwbRq1Ypz587FHevfvz9lypShfPnyzJkzBzDTlx966CEqV67MsGHD4s7t3r0769evT7Jci+Vm0mfJQmCjRjzw/vs8tXAhvc6fp/3ff1O3f39ylizJzokTmd25M9+VKsXwwoX5vVUrVn76KQcXLSLCmYyU2rht0h8ReeKmonIich7YrKrxZ6S5C0TkF6C8s5sLOKeqoSJSAtgO7HSOrVDVHu6W7y2yZs3Kli1buHr1KpkzZ+bPP/+kiMtc/pYtW9KyZUuv6dOkSRP69++Pv78/b7zxBv379+fTTz9l27ZtTJgwga1bt3LkyBEaN27Mrl27mDNnDtWqVWPWrFlUrVqVF198kY0bNxITExOXq8JicQf+GTNS5IEHKPLAA9R6801ioqM5vXUrh/7+m6MrVnB0xQr2/P47AOLnR95KlQioVYuAmjUpVLMm+SpVSvErXidFu2cxuYgXOvsNgBUYg/yBqo51p0Kq+lTsexH5HDjvcniv2wcJF74MJza4tUkKhELDwbc9rXnz5sycOZM2bdrw888/0759e/766y8AxowZw5o1a/jqq6/Yv38/HTp0ICoqioceeiiuvqrSq1cvFixYQMmSJRPM4TB37lz69u1LeHg4pUuXZvTo0bf0qps2bRr3vnbt2nE964TSaqZPn56rV68SFRUVV++9996zs94sHscvXTryBweTPziY+158EYCrp09zdNUqjq5cydEVK9g9aRKbv/sOAP/MmSlYtSqFatakUI0aFKpRg1ylS6eoXBlJmeIcA1RU1daq2hoIAsIxeYrf8JRiYu5SW+BnT8nwNe3atWPChAlcu3aNTZs2xaXCvJnevXvzwgsvsHr1agoVKhRXPmXKFHbu3MnmzZsZOXIky5Ytu6XuqVOn+Oijj5g3bx7r1q2jevXqfPFF4onzRo0aRfPmzQGTb6JYsWJxx2LTajZp0oRjx45Rq1YtXn/9daZNm0a1atUSTDBksXiSzHnzUqp5cx7o1482f/xBz9OneXb3bh7+6SeCu3c3uU+++YaZHTrwfdmyfJUnDxMbN2bRa6+xbdw4Tm7ZQoxLp8LbJKVHXEJVXfMWngDKqeoZEfFkto+6wHFV3e1SVlJE1gMXgHdV9a/4KopId6A7QGBgYOJSktBz9RTBwcGEhYXx888/8/DDDyd43tKlS5k0aRIATz/9NG+8YZ5/S5YsoX379qRLl47ChQvTqNGtufpXrFjBtm3b4hKpR0REJJok/uOPP8bf35+OHTsCCafV9Pf3Z/z48YDJJNesWTOmTZvGK6+8wsGDB+ncubNXXSsWiysiQu4yZchdpgxBznc5OjKS01u3cmzNGo6tXs3xtWtZP3Ro3EBguowZyVe5MvlDQigQEkL+kBDyBweTKXduj+ubFEP8l4jMAH519lsDS0QkK3DuToSKyDygUDyH3lHV2Ewx7bmxN3wUCFTV0yJSDfhdRCqp6oWbG1HVEcAIMFETd6Kjt2jZsiWvvvoqixYt4vTphNdjTehv1O3+XqkqTZo04eefb//H4ocffmDGjBnMnz8/rt2kpNUcNmwYXbp0Yfny5WTIkIFffvmFOnXqWENsSVGkS5+eAqGhFAgNJfi55wCIiYri9I4dnNywgRMbN3Jywwb2Tp/OllGj4uplDwyMc4Xkr1KFfMHB5ClXzq26JcUQ98QY3wcAAX4EJjlrNt3R6h2q2jix487svSeAai51wjEuEVR1rYjsBcoBqTo2rVu3buTMmZMqVaokuELFAw88wIQJE+jUqRPjxo2LK69Xrx7ffvstnTt35sSJEyxcuJAOHTrcULd27dr07NmTPXv2UKZMGa5cucKhQ4cod9MX6Y8//uDTTz9l8eLFcQnfwTwoOnTowCuvvMKRI0duSat59uxZZsyYwdy5c5k2bRp+fn6ICNeuXXPD3bFYPIufv79ZRqpyZYI6dQJM5+XysWOc3LSJkxs3xr2G/fFHnPsinUt+bXeQlKWSFPjN2bxFY2CHqh6KLRCR/MAZVY0WkVJAWWCfF3XyCEWLFqV3796JnjNkyBA6dOjAkCFDaN26dVx5q1atWLBgAVWqVKFcuXLUr3/rHJv8+fMzZswY2rdvH7dixkcffXSLIX7ppZcIDw+nSZMmgDHgw4cPvyGtpr+//y1pNT/44APeffddRIRmzZrx9ddfU6VKFXr0SLUBLZY0joiQLSCAbAEBlGzWLK48KjycMzt2cGrzZk5u2gQDB7pPZhImdDwBfAoUwPSIBWOfk762T3KVEhmDCU8b7lLWGvgAiAKigb6qOv12bdkJHRaLxRO4c0JHUlwTnwGPqqrXVmtU1WfiKZsETPKWDhaLxeItkhK+dtybRthisVjSGknpEa9xZrv9jjNYBqCqkz2llMVisaQlkmKIcwBXgKYuZQpYQ2yxWCxuIClRE129oYjFYrGkVZKS9CcTJt9EJSBTbLmqdvOgXhaLxZJmSMpg3VjMLLhmwGKgKHDrGuuWZHHu3DnatGlDhQoVqFixIsuXLwfgzJkzNGnShLJly9KkSRPOnj0bb/0NGzZQu3ZtQkNDqV69OqtWrYo7Fl/ayptp0KAB5cuXJyQkhAceeICdO3fGe54nCQsLi5smDbBmzRr+85//ACbp0UsvveR1nSwWn6CqiW6Y9ekANjmv6YEFt6uXUrZq1appSqRz5846cuRIVVUNDw/Xs2fPqqrqa6+9pv3791dV1f79++vrr78eb/0mTZrorFmzVFV15syZWr9+fVVV3bp1qwYHB+u1a9d03759WqpUKY2Kirqlfv369XX16tWqqvrtt9/qo48+miS9Y2JiNDo6OsnXmRgLFy7UFi1axHts9OjR2rNnT7fIsVg8AbBG3WSnkjJYF5vY55yIVAaOASXc/kTwFcdfhmsb3NtmplAoODjBwxcuXGDJkiWMGTMGgAwZMpAhQwbApJ2MnercpUsXGjRowKeffnpLGyLChQsmzcb58+fj8j8klLYysUQ/9erVY/Bgo+/AgQOZOHEi4eHhtGrVivfff5+wsDCaN29Ow4YNWb58Ob///ju//PILY8eOxc/Pj+bNmzNgwAD27t1Lz549OXnyJFmyZGHkyJFUqFCBZ555hhw5crBmzRqOHTvGZ599Rps2bXjzzTfZvn07oaGhdOnShfvuu49BgwYxY8aMG/Q7efIkPXr04ODBgwAMHjw4LomRxXIvkBRDPEJEcgPvAdOAbMB/ParVPc6+ffvInz8/Xbt2ZePGjVSrVo0hQ4aQNWtWjh8/TkBAAAABAQGcOBF/7v3BgwfTrFkzXn31VWJiYuJSYB4+fJjatWvHnRebtjIxpk+fTpUqVZg7dy67d+9m1apVqCotW7ZkyZIlBAYGsnPnTkaPHs2wYcOYPXs2v//+OytXriRLliycOXMGMKtzDB8+nLJly7Jy5UpefPFFFixYAMDRo0f5+++/2bFjBy1btqRNmzYMGDDgBsObUK6N3r1706dPHx588EEOHjxIs2bN2L7dhrZb7h2SEjXxnfN2MVDKs+r4gER6rp4iKiqKdevWMXToUGrVqkXv3r0ZMGAAH374YZLb+Oabb/jf//5H69atmThxIs8++yzz5s1LMG1lfHTs2JHMmTNTokQJhg4dypAhQ5g7d27cChuXLl1i9+7dBAYGUrx48TgDP2/ePLp27RqXHChPnjxcunSJZcuW8eSTT8a1H5vbAuDxxx/Hz8+PoKAgjh93zap6e+bNm8e2bdvi9i9cuMDFixfJnj17stqxWFIqSYmayIjJvlbC9XxV/cBzat3bFC1alKJFi8Ylgo/tHQIULFiQo0ePEhAQwNGjRylQoAAAXbt2Zf369RQuXJhZs2bxww8/MGTIEACefPJJnnPS+iUlbWUs48aNo3r161PlVZW33nqL559//obzwsLCyJo16w3n3WzcY2JiyJUrFxs2bIhXVkaXbFXxPSwSIyYmhuXLl5M5c+Zk1bNYUgtJiZqYCjyGSbZz2WWz3CGFChWiWLFicZEK8+fPJygoCDBpJ3/44QfA5Ad+7LHHABg9ejQbNmxg1qxZABQuXJjFixcDsGDBAsqWLRtXf8KECYSHh7N///5b0lYmRrNmzRg1ahSXnAUYDx8+HK9rpGnTpowaNYorV64AJtIjR44clCxZkl9/NWmrVZWNGzcmKi979uxcvHj7AJymTZvy1Vdfxe0nZOwtltRKUnzERVX1odufZkkOQ4cOpWPHjkRERFCqVClGjx4NwJtvvknbtm35/vvvCQwMjDNsNzNy5Eh69+5NVFQUmTJlYsSIEQC3TVuZGE2bNmX79u1xA3vZsmXjp59+uqX+Qw89xIYNG6hevToZMmTg4Ycf5pNPPmHcuHG88MILfPTRR0RGRtKuXTtCQkISlBccHIy/vz8hISE888wzCS46+uWXX9KzZ0+Cg4OJioqiXr16dm08yz1FUtJgjgCGqupm76jkXmwaTIvF4gm8kgZTRDZjckr4A11FZB8m6U9sPuJgdyhgsVgsaZ3EXBOPeE0Li8ViScMkOFinqgdU9QAQgFmiKHb/DPEv/GmxWCyWOyApURPfAJdc9i87ZRaLxWJxA0kxxKIuI3qqGkPSoi0sFovFkgSSYoj3ich/RCS9s/XmHlg92WKxWFIKSTHEPYD7gcPAIaAW0N2TSqUF0qVLR2hoKJUrV+bRRx/l3LlzABw5coQ2bdokWrdEiRKcOnXqruSfO3eOYcOGxXvsn3/+oWHDhlSsWJFKlSrFzeCDhNN0Ll26lODgYGrUqMGePXviZDRr1izZM+ksljSHu9K4pdQtpabBzJo1a9z7zp0760cffZTkusWLF9eTJ0/elfz9+/drpUqV4j125MgRXbt2raqqXrhwQcuWLatbt25V1YTTdLZq1Up37dqlc+fO1VdeeUVVVV955RVdtGjRXelpsaRU8HIazHuaKxc3EB15zq1tpkufiyzZQ5N8fp06ddi0aRNg8jo88sgjbNmyhejoaN544w3mzJmDiPDvf/+bXr16AWZm3vTp04mMjOTXX3+lQoUKXL58mV69erF582aioqLo168fjz32GFu3bqVr165EREQQExPDpEmTeO+999i7dy+hoaE0adKEgQMHxukTEBAQlwEue/bsVKxYkcOHDxMUFJRgms706dNz9epVrly5Qvr06dm7dy+HDx+mfv367rmpFss9TJo3xL4mOjqa+fPn8+yzz95ybMSIEezfv5/169fj7+8fl24SIF++fKxbt45hw4YxaNAgvvvuOz7++GMaNWrEqFGjOHfuHDVr1qRx48YMHz6c3r17x02pjo6OZsCAAWzZsuW2eRvCwsJYv359XIKihNJ0vvXWW3Tv3p3MmTMzduxYXn311WRlk7NY0jJp3hAnp+fqTq5evUpoaChhYWFUq1aNJk2a3HLOvHnz6NGjB/7+5mPKkydP3LEnnngCgGrVqjF5sllQe+7cuUybNo1BgwYBcO3aNQ4ePEidOnX4+OOPOXToEE888URcgqDbcenSJVq3bs3gwYPJkSNHoueGhoayYsUKAJYsWULhwoVRVZ566inSp0/P559/TsGCBZMk12JJayRlsA4AEaktIgtEZKmIPO5BndIEmTNnZsOGDRw4cICIiAi+/vrrW87ReNJNxhKbVjJdunRERUXFnT9p0iQ2bNjAhg0bOHjwIBUrVqRDhw5MmzaNzJkz06xZs7hk7YkRGRlJ69at6dixY5zRh+tpOoEb0nS66vzRRx/x3nvv8f777/P+++/TqVMnvvzyy6TdGIslDZKgIRaRm2fPvQK0BB4C7H9ON5EzZ06+/PJLBg0aRGRk5A3HmjZtyvDhw+MMratrIj6aNWvG0KFD46IU1q9fD5gVQUqVKsV//vMfWrZsyaZNmxJNQamqPPvss1SsWJFXXnnlhmMJpemM5YcffqBFixbkzp2bK1eu4Ofnh5+fX1zKTIvFEg8JjeIBv2OWR8rk7I/AhK09Byx112ihp7fUEDWhqvrII4/ojz/+eEM0Q2RkpPbp00crVqyowcHBOnToUFW9MWpi9erVcQuHXrlyRbt3766VK1fWSpUqxS3M+cknn2hQUJCGhIRos2bN9PTp06qq2r59e61UqZK++uqrN+jy119/KaBVqlTRkJAQDQkJ0ZkzZ6qq6qlTp7RRo0ZapkwZbdSoUVxbqqqXL1/WBg0aaEREhKqqLlmyRCtXrqxVq1bVnTt3uvP2WSw+BzdGTSSaBlNEHgV6Az8Ak4AOQBbgZ1U96bGngxuxaTAtFosncGcazER9xKo6HWgG5AImAztV9cvUYoQtFoslNZCYj7iliPwNLAC2AO2AViLys4iU9paCFovFcq+TWPjaR0AdIDMwS1VrAq+ISFngY4xhtlgsFstdkpghPo8xtpmBuBUkVXU31ghbLBaL20jMR9wKMzAXhRmks1gsFosHSLBHrKqngKFe1MVisVjSJEmeWedORORJEdkqIjEiUv2mY2+JyB4R2SkizVzKq4nIZufYl5LQlDOLxWJJZfjEEGOiMJ4AlrgWikgQxv9cCTODb5iIpHMOf4OZUFLW2R7ymrYWi8XiQXxiiFV1u6rujOfQY8AEVQ1X1f3AHqCmiAQAOVR1uTOj5Ufgce9pbLFYLJ4jpWVfKwKscNk/5JRFOu9vLo8XEenO9VVEwkVki5v1TCr5gLtbSsPKtrKt7JQqu7y7GvKYIRaRecDNiYMA3lHVqQlVi6dMEymPF1UdgcmNgYiscdc0xORiZVvZVva9LdtdbXnMEKtq4zuodggo5rJfFDjilBeNp9xisVhSPb4arEuIaUA7EckoIiUxg3KrVPUocNHJiSxAZyChXrXFYrGkKnwVvtZKRA5hplDPFJE5AKq6FZgIbAP+AHqqarRT7QXgO8wA3l5gdhLFjXCn7snEyrayrWwr+7YkmgbTYrFYLJ4npbkmLBaLJc1hDbHFYrH4Gnct9eGtDRNVsRDYDmwFejvleYA/gd3Oa26nPK9z/iXgq5vaqgZsxvidv8Rx1bhRdhNgrSNjLdDIi7JrAhucbSPQyluyXeoFOvf9VS9edwngqsu1D/fmdQPBwHLn/M1cX2rM09fd0eWaNwAxQKiXZKfHrOKz2anzlhc/7wzAaEfGRqCBB2Q/6ezHANVvqvOW0/5OoNkdy07sYErcgACgqvM+O7ALCAI+A950yt8EPnXeZwUeBHpwqyFehRkwFMzgX3M3y74PKOy8rwwc9qLsLIC/S90TLvsele1SbxLwKzcaYk9fdwlgSwJteVq2P7AJCHH28wLpvHnPnfIqwD4vXncHzIzY2O9dGFDCS7J7AqOd9wUwHR4/N8uuiJm8sQgXQ+wc2whkBEpiggju7PNO7GBq2DBhbE0wT6QAlxu686bznsHFEDvn7HDZbw986wnZTrkAp50PzduySwLHMYbCK7IxU9AHAv1wDLE3ZJOAIfaS7IeBn3wh+6ZzPwE+9uJ1twemO9+vvBgDlsdLsr8GOrmcPx/zb9Btsl32F3GjIX6LG3v/czDGN9myU7WPWERKYHqdK4GCauKNcV4L3KZ6EZIxbdoNslsD61U13FuyRaSWiMT+Re6hqlHekC0iWYE3gPdvqu6te15SRNaLyGIRqetF2eUAFZE5IrJORF73omxXngJ+9qLs34DLwFHgIDBIVc94SfZG4DER8XfmHlTDuBjcKTshigD/xCMj2bJTWq6JJCMi2TB/fV9W1Qt3kBUzWdOm70a2iFQCPgWaelO2qq4EKolIReAHEZntJdnvA/9T1Us3neMN2UeBQFU9LSLVgN+d++8N2f4YN1gN4AowX0TWAhe8IDv2/FrAFVWNza/ijeuuCUQDhYHcwF9OigNvyB6FcR2sAQ4AyzCLWbhNdmKnJiAj2bJTZY9YRNJjbtQ4VZ3sFB93srThvJ5IqL7DHU2bTq5sESkKTAE6q+peb8qORVW3Y3oslb0kuxbwmYiEAS8Db4vIS96QrSZz32nn/VqM366cl677ELBYVU+p6hVgFlDVS7Jjacf13nCsTp6W3QH4Q1UjVfUEsBSo7g3Zqhqlqn1UNVRVH8OsOL/bzbITwm0pGVKdIXamOH8PbFfVL1wOTQO6OO+7cJsp0HoH06aTK1tEcgEzMX6kpV6WXVJE/J33xTGDDWHekK2qdVW1hKqWAAYDn6jqV1667vzi5LAWkVKYafL7vCEb4yMMFpEszr2vD2zzkmxExA8zwj8htsxLsg8CjcSQFaiN8ZF64/PO4shERJoAUarq7nueEO5LyZAc53VK2DB//RQzOr3B2R7GDBLMxzwN5wN5XOqEAWcwoVSHgCCnvDomSf1e4CtuH96SLNnAu5ie6AaXrYCXZD+NCbnZAKwDHndpy6Oyb6rbjxujJjx93a2d697oXPej3rxuoJMjfwvwmZdlNwBWxNOWp+95Nkx0zFZMeoLXvCi7BGYgbzswDyjuAdmtMHYjHDPoPcelzjtO+ztxiYxIrmw7xdlisVh8TKpzTVgsFsu9hjXEFovF4mOsIbZYLBYfYw2xxWKx+BhriC0Wi8XHWENsuedx4lv/FpHmLmVtReQPX+plscRiw9csaQIRqYyJdb0PSIeJEX1Ir892TE5b6fT6El4Wy11jDbElzSAin2Em2GR1XotjUkb6A/1UdaqYZC9jnXMAXlLVZSLSAOiLyWURqqpB3tXeci9jDbElzeBMhV0HRAAzgK2q+pMzFX0VpresQIyqXhORssDPqlrdMcQzgcqqut8X+lvuXVJt9jWLJbmo6mUR+QUz1b0t8KiIvOoczoRZUeQI8JWIhGIyipVzaWKVNcIWT2ANsSWtEeNsArRW1Z2uB0WkHyafQAhmMPuay+HLXtLRksawUROWtMocoJeTHQsRuc8pzwkcVdUYTOKkdD7Sz5KGsIbYklb5ELPo5SYR2eLsAwwDuojICoxbwvaCLR7HDtZZLBaLj7E9YovFYvEx1hBbLBaLj7GG2GKxWHyMNcQWi8XiY6whtlgsFh9jDbHFYrH4GGuILRaLxcf8P8PEDIG4m5U+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot results\n",
    "years = df1['year']\n",
    "\n",
    "colors = ['darkred', 'red', 'darkorange', 'gold', 'palegoldenrod']\n",
    "\n",
    "for i, q in enumerate(qlabs):\n",
    "    plt.plot(years, d[q], c=colors[i])\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('% change in GDP pc relative to reference')\n",
    "plt.title(f\"GDP pc in {SSPtest} relative to {SSPref} (avg across IAM)\")\n",
    "plt.legend(qlabs, frameon=False, loc=\"lower left\")\n",
    "plt.hlines(0, min(df1['year']), max(df1['year']), colors='black', linewidth=.6)\n",
    "plt.ylim([-100, 25])\n",
    "plt.locator_params(axis=\"y\", nbins=5)\n",
    "plt.rcParams[\"figure.figsize\"] = (5,4)\n",
    "plt.margins(0)\n",
    "\n",
    "if savefig:\n",
    "    if not os.path.isdir(output_dir):\n",
    "        os.makedirs(output_dir)\n",
    "    plt.savefig(os.path.join(output_dir, '{}_income_compare_{}.pdf'.format(SSPtest, SSPref)))\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24543ad1-edbb-48c6-9c7b-b67c7cb87089",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
